% This file is used to define the Global variables, initial conditions and call the required files.
% The variables are:
% p: represents zeta : wave height
% q : represents phi : surface velocity potential
% 1-D case .. in X direction


%--- main parameters
global U  h  g  p  q  nt dx h1 x tmax L steps y itermax Q p_global xrange dx1 dx2;
range=[ 1 2 4];
for k = range
    %clear all;
    
    %----- Physical Parameters----------
    
    U = 1; % Current Velocity of the wave. Assumed Constant
    h = 50; % height of the bed, assumed constant
    g=  9.81; % Gravitational Constant
    
    %index = find(range==k);
    %------------ Simulation Parameters % -----------------
    tmax = 100;     % max time, s
    L=  100;
%     
%     dx_a = 1.0;
%     dx_b = 0.5;
%     x1=linspace(0,L/4,26);
%     x2=linspace(L/4+dx_a,3*L/4,100);
%     x3=linspace(3*L/4+dx_b,L,25);
%     
    d2 = L/10;
    d1 = (L-d2)*0.5;
    d11 = d1+d2;
    d12 = d11 + d1;
    
    N1 = 46;
    dx_a =d1/(N1-1);
    x1=linspace(0,d1,N1);
    N2=k*10;
    dx_b = d2/N2;
    
    x2=linspace(d1+dx_a,d11,N2);
    x3=linspace(d11+dx_b,d12,N1-1);
    
    
    xrange= [x1 x2 x3];
%    xrange = linspace(0, L , 101);
    
    
    x = length(xrange);
    
    steps = x-1;    % number of axial steps
    % Check the CFL condition:
    dx = zeros(steps,1);
    for i = 1:steps
        dx (i) = xrange(i+1)-xrange(i);
        dx1(i) =1/dx(i);
        dx2(i)=dx1(i)*dx1(i);
    end
    
    dt1 = 0.5*min(dx)/sqrt(g*h);
    %--- Initial Conditions ----------
    pi1=1;
    qi=1;
    
    for i=1:x
        p(i,1) = pi1*cos(4*pi*xrange(i)/L);
%        p(i,1) = 0.25*pi1*cos(4*pi*(i-1)*dx/L) + 0.25*pi1 *cos(8*pi*(i-1)*dx/L)  + 0.25*pi1*cos(16*pi*(i-1)*dx/L)+0.25*pi1*cos(32*pi*(i-1)*dx/L)  ;
        q(i,1) = qi;
    end
    
    % Run the matlab based Codes
    i=1;
    
    Complete_Matlab(2);
    A(1:x,i)=y(end,1:x);
    B(1:x,i)=y(end,x+1:2*x);
    % i=2;
    % Complete_Matlab(i);
    % A(1:x,i)=y(end,1:x);
    % B(1:x,i)=y(end,x+1:2*x);
    
    % Run the custom Program
  % itermax=9;
    %dt =0.023;
   % CFL = [1.8 0.9 0.4];
    %CFL=1.9;
    %dt_v = CFL*0.5*dx/sqrt(g*h);
    %dt_v=[CFL*dt1];
    %dt_v = dx*50;
    dt_v = [0.01];
    for i = 1:length(dt_v)
        dt = dt_v(i);
    %    X = sprintf(' Simulation with Nodes: %f , CFL %f ,  dt %f',x,CFL(i),dt);
    %    disp(X);
        TimeDiscretizers(dt,1);
        A(1:x,i+1) = p;
        B(1:x,i+1)=q;
        X = sprintf('The Norm of the Error for Zeta is %f' , norm_L2(A(:,1),A(:,i+1)));
        disp(X);
        %X = sprintf('The Norm of the Error for varphi is %f' ,  norm_L2(B(:,1),B(:,i+1)));
        %disp(X);
    end
    plot(xrange,A(:,1:2));
    pause
   %     A_G(:,:,index) = A;    
end


